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Abstract. The spatio-temporal dynamics of a population present one of the most fascinating as- 
pects and challenges for ecological modelling. In this article we review some simple mathematical 
models, based on one dimensional reaction-diffusion-advection equations, for the growth of a pop- 
ulation on a heterogeneous habitat. Considering a number of models of increasing complexity we 
investigate the often contrary roles of advection and diffusion for the persistence of the population. 
When it is possible we demonstrate basic mathematical techniques and give the critical conditions 
providing the survival of a population in simple systems and in more complex resource-consumer 
models which describe the dynamics of phytoplankton in a water column. 
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1. Introduction 

For long times field biologists and naturalists have been fascinated by the richness and beauty 
of complex patterns that can be observed in spatially extended populations. However, the same 
observations also constitute a challenge for theoreticians who aim to explain this complexity by 
means of mathematical models. At first glance one might be tempted to argue that the spatial 
diversity of natural populations mainly originates from some underlying abiotic heterogeneity of 
the environment. If growth conditions vary between different locations then this spatial variation 
should be reflected in the density distribution of natural populations. Thus, it is reasonable to 
assume that a large part of the observed richness in the patterning of biotic landscape can be 
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attributed to the spatial heterogeneity of growth conditions. On the other hand in nearly all systems 
populations at different locations are coupled and are able to interact with each other, be it by 
random, diffusive mixing or by directed, advective flows or currents, which are able to transport 
individual organisms from one region to another. Note that the term "advection" here is used in a 
very general sense, including different mechanisms of transport, such as drift caused by a current 
of water, sinking and floating-up in the gravity field, chemotaxis, etc. 

The impact of such spatial interactions on the fate of a local population can be very diverse, 
depending on whether the flow brings-in new immigrants into the habitat or if it takes them away. 
Furthermore, diffusion and advection can have opposing influences, where, for example, one pro- 
cess may be beneficial for the population, while another process has adverse effects. Usually, the 
role of an advective flux is evident and mainly depends on its strength and direction. In contrast, 
the role of diffusivity can be of two kinds. On the one hand, diffusivity accelerates the spread of 
a population in a habitat and is necessary to provide the population's persistence in a flux. On the 
other hand, intensive mixing may transport too many organisms into unfavourable zones, result- 
ing in the extinction of the population. Moreover, in consumer-resource models usually also the 
resource fluxes are driven by diffusivity and advection, a fact which leads to even new patterns 
and dynamical behaviour. As a consequence, it is quite possible to find stable populations in loca- 
tions where growth conditions alone would not permit persistence. On the other hand, seemingly 
well-being habitats may not be able to sustain a stable population, if they suffer population out- 
flows. From all these effects, populations are spatially structured in an intricate interplay between 
local growth and its geographical variations on the one hand and spatial transport by diffusion and 
advection on the other hand. 

This article is devoted to present an introductory review of these topics and to introduce the 
reader some of the most commonly used models for the population dynamics in a non-uniform 
turbulent environment. Even though we have in particular the specific case of phytoplankton dy- 
namics in mind, the main concepts considered here hold for any population subjected to mixing 
and advection. When it is possible, to reproduce the whole picture we perform derivations of the 
main mathematical relations. Otherwise we aim at least to demonstrate some basic ideas and refer 
to detailed discussions, providing information about the main techniques that are useful in this 
field of research. There exist many excellent reviews about spatial population dynamics (see e.g. 
Il32l |47l [651 1721 ). Nevertheless to our knowledge, the interplay of advection and diffusion for a 
heterogeneous population, as elaborated in this text, has never been described. So, while the text 
will not provide much new insights for the specialist, we hope that many readers will find this text 
a valuable introduction into the basic mechanism and critical conditions providing the survival of 
a population in a heterogeneous environment. 

To focus on the main ideas, we had to confine the review to topics related with population 
survival. However, we would like to briefly mention other important aspects of spatial population 
dynamics that had to be omitted here. First, we had to restrict to single species populations, while 
one of the most important and still debatable challenges concerns the competition of spatially 
structured populations H104[|25ll98l and the high diversity of phytoplankton species (see e.g. If88l 
and the references therein). Second, for the sake of simplicity we included only Fickian diffusion 
models, whereas in the ocean diffusivity depends on the scale of phenomena [1721 . Third, we had 
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to omit the interesting aspects concerning conditional persistence in situations where the growth 
is influenced by an Allee effect [1 10011 . Fourth, the interaction of turbulent currents with other pro- 
cesses can produce complex spatial structures in phytoplankton distributions ||57ll , whereas locally 
oscillatory behaviour may lead to spatio-temporal chaotic waves Ir79ll80l . Finally, we did not even 
touch such important problems as model validation and the simulation of multicomponent natural 
ecosystems 11601 . 

The review is structured as follows. In chapter two, we shortly discuss some non-spatial mod- 
els, which will be used later to build-up spatial explicit models. In chapter three we consider 
critical patch models, starting from the well known KiSS model, and proceed with more sophis- 
ticate models, which draw out the influence of boundary conditions, finiteness of the favourable 
patch, heterogeneity of the mixing, etc. Furthermore, we consider the Fisher- Kolmogorov equation 
and its extensions to advection, non-uniform mixing and heterogeneous environments. In the last 
chapter, we consider ID models of the vertical phytoplankton distribution, starting from those that 
include only light limitation and then considering models including many limiting factors. How- 
ever, even in the latter models we focus on the basic theoretical aspects and general conclusions. 

2. Non-spatial models 

Logistic growth As a warming-up, in this chapter we shortly describe some simple population 
models in a non-spatial context. These models will mainly be used as a base for developing spa- 
tially extended models, but they may also serve as a benchmark for exploring the properties of their 
spatial analogues. Let P{t) denote the density of a population of interest at time t. In the simplest 
way, the dynamics of P can be modelled in terms of an ordinary differential equation of first order 



where the growth rate /u(P) depends on the population density P. The function /i(-P) takes into 
account for density dependent regulatory factors which usually are associated with resource deple- 
tion or conditions of over-crowding. Consequently, (neglecting the possibility of an Allee effect) 
we may assume that fi(P) is monotonically decreasing and eventually becomes negative for large 
densities. The simplest form to put this into a model is the logistic growth 



where K is the carrying capacity of the system and describes the maximal population density that 
can be sustained by the system, i.e. n(K) = 0. Equation (12.21) can be thought to arise from a 
Taylor expansion of equation (12.11) for small densities P <C K, where /io = /i(0) is the growth rate 
at small density and K = —fi (d/i/dP^ 1 is the carrying capacity. 

Resource-limited growth As a secondary model class, we investigate resource limited popula- 
tion growth, where the dynamics of the resource is explicitly included into the model equations. 



dP 
~dt 



(2.1) 




(2.2) 
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Figure 1 : Typical phase trajectories of the non-spatial models. (A) the dynamics of the conservative 
model (12.31) is reduced to a straight line; (B) the non-conservative model (12.71) leads to the extinction 
of the population; (C) the trajectory of system (12.101) spirals in phase space to an isolated fixed 
point. The initial values P(0) = 0.1, N(0) =1.5 are marked by filled circles, other parameters are 

c = 5, H N = 5, m = 1, S = 0.005, e = 0.5, and Ni = 15. 



Possibly the most simple way to model an isolated consumer-resource system (i.e., in the absence 
of any fluxes out of the system) is given in the following form 

dP 

— = fi(N)P-mP, 

(2.3) 

— = -u(N)P + mP . 
at 

Here N(i) denotes the resource concentration (e.g., a limiting nutrient such as nitrate or phosphate) 
at time t and the population abundance P(t) is measured in units of the organism's resource con- 
tent. As usual, in contrast to equation (12.11) . here we explicitly separated the growth into a term 
proportional to the resource uptake, fi(N)P, and a death process with mortality m. Note that in this 
review the symbol fi has a somewhat different interpretation in the context of a resource-implicit 
(12.11 ) or resource-explicit model (12.3b . We hope that the reader is not confused by this notation. 

Usually, [x(N) is assumed to vanish for small resource concentrations (p(0) = 0) and to be 
saturated for large N. A convenient parametrisation is given by the Monod form 

with the half-saturation constant Hjq. 

Equation (12.31) represents a closed system, where recycling is assumed to be 100% efficient. 
Thus the biomass that is taken out of the system due to death processes is remineralised into the 
nutrient pool, i.e. no material is lost in the conversion from consumer to nutrients. As a conse- 
quence, this system can be reduced to a one-dimensional model, despite the fact that it contains 
two equations. The reason is that the total resource concentration, either free or bound to the pop- 
ulation, is a first integral of motion: S = N + P = const. After plugging this into the equation for 
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P we obtain a single equation of the form (12.11 ) for the time dynamics of the population density, 

dP 

— = P fr(S -P)-m). (2.5) 

If the saturation of the uptake rate is neglected, /x(iV) = cN, we retrieve the logistic growth model 
(TO) with 

jiQ = cS — m and K = S . (2.6) 

c 

The phase trajectory of system (12.31) in (P, N) coordinate is a straight line from the initial condi- 
tions to an equilibrium state (Fig. UK). 

We note though, that the reduction to a single equation works only for a closed system. In 
general, recycling efficiency will never be fully effective. Thus, one may assume that only a 
fraction e of the dead biomass is remineralised, leading to a non-conservative model 

^ = fi(N)P - mP , 

dt (2.7) 
= -n(N)P + emP . 

The two approaches introduced in this section will be used in chapter 4 as building blocks to 
model the local reaction of a spatial extended population. However we should note that both 
(12.31 ) and (12.71 ) have some shortcomings, which ultimately are related to the fact that the models 
do not contain in- and out-flows. For example, the model (12.71) does not even admit isolated 
stationary states in the (N, P)-phase plane and for any initial condition leads to the extinction of 
the population, P(t) — * (Fig.QJS). While this is no problem in the spatial extended counterparts 
of the models, which necessarily have exchanges with surroundings specified by the boundary 
conditions, in local models these exchanges should be introduced as additional terms. 



Chemostat models To show that these problem disappear as soon as external fluxes are incorpo- 
rated into the models, we shortly introduce the most simple and frequently used chemostat model. 
Consider a well- stirred reactor that contains phytoplankton cells with concentration P(t) and a 
limiting nutrient with concentration N(t). The chemostat is supplied with the nutrient at an input 
concentration iVj from an external nutritive medium. The outflow contains both medium and phy- 
toplankton cells. Inflow and outflow are characterised by the dilution rate 5. Then the chemostat 
equations take the following form 

/j(N)P - (m + 5)P , 

(2.8) 

S(Ni — N) — fi(N)P + emP . 

In this model, the total amount of free and bound nutrients S = N + P is not conserved. Thus, the 
chemostat sustains a stable feasible equilibrium and the system spirals in phase space to an isolated 



dP 
~dT 
(W 
~dT 
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fixed point (Fig.[T£), see e.g., ll38l[T0ll39l . Furthermore, for perfect nutrient recycling (s = 1) one 
can easily show that S(t) follows the simple dynamics 

^ = 5(iV, - S) . (2.9) 

Therefore, after some transient time, determined by the exchange rate 5, the system settles to a state 
where S = iVj. This means that the total amount of nutrients in the system, either bound to the pop- 
ulation or free, asymptotically goes to the external nutrient concentration, and one asymptotically 
retains a one-dimensional equation, similar to (12.51) 

dP 

— = P[ f i(N l -P)-(m + 5)]. (2.10) 
at 

3. Critical patch models for an extended population 

In this section we review some simple one-dimensional spatio-temporal models which demonstrate 
the intricate interplay between mixing, advection, boundary conditions and other factors which are 
critical for the persistence of a population. Consider a population growing on a one-dimensional 
(ID) axis, which may represent either a horizontal or vertical spatial extension. Let P(x, t) denote 
the population density at time t and position x. We assume that the change of density occurs as a 
result of the local death-birth processes and of the spatial movement of organisms due to diffusion 
and advection. The dynamics of such a population can be written in terms of a reaction-diffusion- 
advection equation [|72j [621 

dP(x, t) 

= reproduction — advection + mixing 



Of 



dP d dP 

(3-D 

ox ox ox 



where fi represents the growth rate (compare to equation (12.11) ). v is the advection velocity and D 
is the diffusivity, which can depend on x and other variables. The advective term describes 
the drift of the population with the constant velocity v. There are many physical or biological 
factors which can give rise to advective processes. For instance, in a vertical system the organisms 
(e.g. phytoplankton cells) may sink (or rise) because they are more heavy (or light) than their 
surrounding medium. Another example arises for a population in a vertical/horizontal stream (e.g., 
in a horizontal flow of a river or in an ocean current). Advection and diffusion can be written in 
terms of a derivative -§^J{P, x) so that the whole flux of biomass is given as 

dP 

J = vP- D— . (3.2) 

ox 

Additionally, model (13.11) has to be complemented by appropriate boundary conditions which 
specify the environmental influence on the spatially extended system. One can distinguish between 
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three important types of boundary conditions. The first type (or Dirichlet) boundary condition 
merely specifies the value of the solution at the boundary. For example, the condition 



P(0) = 0, P{L) = 



(3.3) 



assumes that the population density vanishes at the habitat edges, which corresponds to an abso- 
lutely hostile environment outside the segment [0, L]. The second type (or Neumann) boundary 
condition specifies the flux across the boundary. For instance, the condition 



J(x) 



x=0,L 



vP -D 



dP 

dx 







(3.4) 



=0,L 



specifies impenetrable boundaries, in other words there is no in- and outflow of biomass across 
the habitat edges. The boundary condition of the third (Robin) type is a linear combination of 
the Dirichlet and Neumann conditions. Consider, for instance, a penetrable barrier, where the flux 
of organisms across the barrier is proportional to the difference of the population density outside, 
Pout, and inside, P(x), yielding 



J{x% 



=0,L 



VP 



^ dx J 



h(P(x) - P t 



out J 



=0,L 



(3.5) 



x=0,L 



If the permeability of the barrier h — this condition is equivalent to the zero flux boundary 
condition, and as h —> oo it approaches the Dirichlet boundary conditions P(L) = P out . 



3.1. Persistence on a finite patch: critical patch size and mixing 

For simplicity, we begin our investigation with models which implement only diffusion but no 
advection. We concentrate on the population dynamics on a finite favourable patch and show that 
in this case diffusion plays a negative role for the populations' persistence. 
In the absence of advection, v — 0, equation (13.11) reduces to 

complemented by appropriate boundary conditions. Perhaps the most important biological ques- 
tion that one may ask in this model concerns the conditions under which the population will be 
able to persist or die out. This question can be trivially answered for an infinite homogeneous 
environment, where the growth rate is independent of the spatial coordinate, fj,(x, P) = fJ>(P)- In 
this case (and assuming that there is no multi stability, i.e. the equation fJ,(P) = has a unique 
positive solution P*) the fate of the population entirely depends on the sign of the linearisation of 
the growth rate around P = 0. The population can invade if and only if the linearised growth rate 
Ho = n{P = 0) > 0. Furthermore, it is easy to verify that the final distribution must be uniform, 
P(x) = P*. 

In reality the growth conditions will differ between various locations so that P) will de- 
pend on x. We refer to such a situation as a "heterogeneous environment". In a heterogeneous sys- 
tem, obviously also the linearised growth rate will be a function of the spatial position, fi = Ho(x), 
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which makes the problem of the population persistence more complicated. Even if locally the pop- 
ulation faces bad environmental conditions, fio(x) < 0, due to the population inflow from other 
locations it may still be able to persist. On the other hand, a favourable patch with fio(x) > is no 
longer save, because diffusion out off the patch can cause additional losses which may even lead 
to local extinction. 



KiSS model One simple and elegant model to study this problem has been independently in- 
troduced by Kierstead and Slobodkin Il44l and Skellam ||96ll . The main idea is to separate the 
landscape into a favourable area (which will be denoted as the species' habitat) adjoining some 
hostile environment from both sides. For its simplicity, Akiro Okubo assigned to the model the 
name KiSS, which on the one hand includes the authors' initials, and from the other hand can be 
deciphered as "Keep it Simple Stupid" ||9~71l . The KiSS model suggests a population growing on 
a finite patch of size L, surrounded by an absolutely hostile environment with infinite mortality 
(Fig. |2]). Thus, per definition the population density outside of the favourable patch equals zero, 
and it is sufficient to consider the dynamics inside the segment < x < L, upon the condition that 
the solution vanishes at the boundaries. To simplify the situation even further, the growth term in 
equation (13.61) is linearised for small density and is set constant /i within the habitat, which yields 
the following equations 

dP d 2 P 

— = fj^P + D—, 0<x<L (3.7) 
at ox* 

P(0) = P(L) = . (3.8) 

This model could represent, for instance, the vegetation patterns of coastal plants growing on some 
island in the ocean, where diffusion of plants takes place via seed dispersal. Seeds which are 
transported out of the island into the water cannot survive, which is equivalent to infinite mortality 
outside of the favourable patch. 

A detailed analysis of the KiSS model can be found in 1196114411721 . Here we just present one 
simple solution. Using the method of separation of variables we assume the existence of a solution 
of the form 

P(x,t) = X(x)T(t) . (3.9) 
Substituting this expression into (13.71) we obtain the eigenvalue problem 

r n {t) = x n T n (t) 

DX'^x) + (/i - X n )X n (x) = (3.10) 

where the eigenfunctions X n (x) should satisfy the boundary conditions (13.81) . i.e. X n (0) = 
X n (L) = 0. The solution to the last equation of system ( 13.101 ) is well known 



Uiq — X n //i — A n 

— — — x + B n cos d — — — x . 
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Figure 2: Influence of diffusion on the persistence of an isolated population in the KiSS model 
(13.71) . (A) D = 3.2 < D$, the population growth is unbounded; (B) D = 3.3 > D , the population 
goes extinct. Parameters are yU = 2, L = 4, K = 400, thus D = 3.24. The horizontal thick black 
line shows the location of the favourable patch. 

It satisfies the boundary conditions only if B n = and the argument of the sine is divisible by Tin 
when x = L. From the last condition we obtain the eigenvalues 

<K 2 n 2 D 

A„. = /i j?— , n = l,2,... (3.11) 

Solving the first equation of system (13.101) . we find that T n (t) = Ce Xnt . The problem is linear, 
therefore using the superposition principle, we can write P(x, t) as the sum 

oo 

P{x,t) = X^ sin (^) e Kt , 

n=l 

where the coefficients A n depend on the initial conditions. 

Thus, with the separation ansatz (13.91 ), we were able to present the dynamics of the full model 
as a sum of simple "modes". The dynamics of mode n is determined by its eigenvalue A n . A 
mode grows if A„ > 0, decays if A n < 0, and can also oscillate if the imaginary part of A„ is not 
vanishing. Moreover, as t — ► oo the full solution approaches to the mode that corresponds to the 
largest (principal) eigenvalue A*, because this mode has the fastest time exponent T{t) ~ e A *. 
From (13.1 II) we find that the largest eigenvalue A* = A x . Therefore, 

P(x,t)— A lS in (f)e^. 

The population grows, and so is able to persist, if Ai is positive. By contrast, for Ai < the 
population dies out exponentially on the whole habitat. Thus using expression (13.1 II) for n = 1, we 
obtain the fundamental relation between the diffusivity, the growth rate and the critical (minimal) 
size of the favourable patch, which provides the survival of the population 




(3.12) 
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Okubo iTTTl |72| demonstrated that this equation can be obtained by means of simple dimen- 
sional analysis. Suppose that L = f(D,fx ). Dimensionally, [D] = m 2 / 's, = 1/s, and 
[L] = m. Thereby, the simplest combination of parameters yields 

L = c 

where c is a non-dimensional constant which equals tt in the ID model. The same expression holds 
in 2D systems Il96l l44l . and it can be shown [|32l |9l that the constant c depends on the principal 
eigenvalue, thereby it represents the geometry of the model. 
Equation (13.121) can also be written in the form 

D<D = ^, (3.13) 

yielding the critical, maximum mixing D under which the population can survive. Fig. [2K and 
Fig. 03 demonstrate the population dynamics for different diffusivities. The population density 
grows exponentially, if the diffusivity is smaller than D (Fig. [2]\) and decays exponentially oth- 
erwise (Fig. [2j3). 

These results highlight the negative aspects of diffusivity for a population and reveal an impor- 
tant ecological insight, namely that of a critical patch size. A finite population under the influence 
of random mixing must be larger than a minimal extension L in order to sustain a stable popu- 
lation. This critical size L simply scales as the square root of the strength of mixing divided by 
the growth rate, ( 13.121) . Despite the model's simplicity the fundamental results ( 13.121 ) and (13.131 ) 
have an important ecological message that prevails in more realistic settings. If a favourable patch 
adjoins unfavourable areas, then the higher the diffusivity, the higher will be the loss rate across 
the boundaries, thus the larger should be the internal area and the growth rate on the habitat to 
compensate for these losses. Moreover, one large patch is better for the persistence of a species 
than two smaller patches of the same total size, since two patches would have four ends, that would 
double the loss rate. Diamond and May [fT5l . McMurtrie lf58ll . Cantrell and Cosner [JSj applied this 
concept of a critical patch size to the design of national parks and natural reserves of optimal size 
and form. 




Logistic growth As the growth of biomass in the KiSS model is not limited, one natural exten- 
sion of the KiSS model is to consider a logistic growth function 

3i>(M)_.. „/, P\. D *P i for <x<L (3.14) 



HP (l - £) 



dt \ K J dx 2 

and the same (hostile) boundary conditions (13.81) . For this model only approximate solutions are 
available (see e.g. Il52l . GO, ll6Tl . Il72~ll ). The stationary solution can be expressed in terms of 
elliptic functions and was investigated by Skellam [|96l , Levandowsky and White [54], and Ludwig 
et al. ||56l . However, applying the invasibility criteria, we can conclude that this system possesses 
the same critical values as the KiSS model. Indeed, as the population density approaches zero, 
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Figure 3: Dynamics of the population with logistic growth, Eq. (I3.14I) . under different values of 
diffusion. (A) Conditions are nearly critical, D = 3.2 < D . (B) Conditions are sufficiently far 
from the survival-extinction transition D = 0.1 <C D . Other parameters are /i = 2, L = 4, 
K = 400, thus D = 3.24. The horizontal thick black line shows the location of the favourable 
patch. 

equation (13.141) goes over to equation (13.71) . This means that if a population can invade in the KiSS 
model, it will also invade in the model with logistic growth, and vice versa. 

Figs. [3] show examples of the population dynamics. In both figures, the resulting population 
density has a characteristic shape, with large densities in the central part of the patch and a decay 
of population numbers, the closer one comes to the hostile border. Note, that if the conditions are 
close to the critical values (13.121) or (13.131) . the growth of biomass becomes limited by the diffusive 
transport, and the maximum of density reached by the population can be much smaller than the 
carrying capacity. This is illustrated in Fig. [3K where max(P(:r)) ~ 6 even though K = 400. 
However, if the conditions are sufficiently far from the critical region, the carrying capacity is 
almost reached in the middle of the patch (Fig. 03). 

To obtain an intuitive understanding into the influence of the density dependence in equation 
(13.141) consider a KiSS model with an effective growth rate of ji* = (yU (l — P(x) / K)) x , which 
equals the average growth rate of the logistically growing population. The new effect now is that 
the effective growth rate ji* decays with an increase of the population density P(x). Therefore, 
also the maximal diffusivity Dq, admitting a survival of the population (13.131) . will be a decaying 
function of P(x). This means, however, that the biomass can grow only until the critical diffusivity 
is reached -Dq(P) = D. Thereby, if initially (when P is negligible) D is close to the critical 
diffusivity of the KiSS model, a small increase of density is enough to decrease and to reach 
the balance between production and loss. 

Finite mortality and other extensions A second unrealistic feature of the KiSS model is the 
assumption of infinite mortality outside of the favourable habitat. This assumption does not hold 
in many important cases and limits the applicability of the KiSS model. For example, it hardly 
could be applied for phytoplankton simulations. Ludwig et al. [56] investigated an extension of 
the KiSS model, where this assumption of an absolutely hostile environment was relaxed. In this 



52 



A. B. Ryabov and B. Blasius 



The role of diffusion and advection 



study the authors assumed a positive growth rate /i inside and a finite mortality m outside of the 
habitat 

{/i , for < x < L 
(3.15) 
— m, for x < or x > L . 

As the mortality is finite, the organisms can survive outside of the habitat and diffuse back. 
This effect reduces the losses at the edges and the population can survive on a smaller favourable 
patch compared to the KiSS model. In the model (13.151) the population can persist if 




(3.16) 



As the mortality approaches infinity this value approaches the critical patch size L of the KiSS 
model. 

Many other examples of critical patch models can be found in the books by Okubo and Levin 
[|72l and Murray [|62l . We just briefly mention some extensions. Okubo ll68l 1701 considered a 
model for growth and diffusion under an attractive force toward the centre of a patch. He found 
that the modified critical size equals L c = L f [v 2 /AaD), where f(0) = 1 and the function 
f(v 2 /AaD) monotonically decreases with v toward zero. Gurney and Nisbet ||27l considered a 
model in which the growth rate parabolically depends on the distance from the habitat centre. This 
approach is more realistic since it includes a gradual transition from favourable to unfavourable 
areas. Wroblewski et al. J 1 1 611 . Wroblewski and O'Brien HI 141 and Piatt and Denman [81] included 
the effect of grazing and obtained an expression similar to (13.121 ) for the critical patch size, in 
which, however, /i is replaced by /i — g, where g characterises the grazing rate. 



Influence of boundary conditions and spatial arrangement It is interesting to extend the anal- 
ysis for more complicated spatial geometries. Seno 11931 and Cantrell and Cosner investigated 
the influence of a spatial sequence of favourable and unfavourable habitats and boundary condi- 
tions. Seno considered a population of organisms migrating between n patches of different quality. 
Cantrell and Cosner [7] used a mean field representation of this problem. They compared the total 
size of a population living on a finite habitat of size L, surrounded by either completely hostile 
environments or by impenetrable boundaries. Using the logistic model (13.141 ), they assumed that 
the habitat possesses patches of different quality, that is, the growth rate changes between positive 
or negative values, fi(x) = ±1, provided that the favourable and unfavourable patches have the 
same total area (Fig. |4]). 

In this model the best spatial configuration of favourable and unfavourable patches depends 
on the boundary conditions. If the exterior region is completely hostile then the location of a 
favourable patch in the middle of the habitat (Fig. UK) provides the best conditions for the popula- 
tion, as two buffer zones separate the favourable region from the completely hostile regions, which 
decreases the population losses. Therefore, the total biomass will be higher than that on the habitat 
shown in Fig. @}3, which in turn is more preferable than the habitat shown in Fig. |4£, where two 
favourable patches adjoin the hostile surroundings. 
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Figure 4: Comparison of several spatial arrangements of favourable (ji(x) = 1) and unfavourable 
patches (fJ,(x) = — 1), assuming either lethal boundaries (top) or impenetrable boundaries (bottom) 
at the borders of the habitat (x — and x — 1) according to 0. Plotted is the final population 
density P(x) in model (13.61 ) assuming logistic growth, fi(x, P) = yu(x)(l — P/4) (solid line). The 
patch quality, fi(x), is shown as dashed line. Each case is classified either as "most favourable", 
"intermediate" or "less favourable", reflecting the total amount of biomass. 



If however, the boundaries act as a barrier, the maximum biomass is achieved if a favourable 
patch directly adjoins one of the impenetrable boundaries (Fig. HJ)). Even though the derivation 
of this result requires rather sophisticated calculations, it is evident that in this configuration there 
is a single favourable patch which has only one border with an unfavourable environment, and 
so all possible losses are minimised. By contrast, any splitting of the favourable patch leads to a 
worsening of the habitat quality (Fig. [4^ and HP). 

It is clear, that the choice of the best arrangement of habitat quality is very important for the 
design of national parks, natural reserves, etc. Further, the same mechanism also captures some 
aspects of relevance for the vertical distribution and competition of phytoplankton species. In a 
water column the surface acts as an impenetrable barrier. Thus, in an incompletely mixed water 
column, a more light limited or buoyant species obtains a competitive advantage to another species, 
whose favourable patch is located in subsurface layers [|34l[89l . 

Arbitrary spatial dependence of the growth rate The behaviour of model (13.61) in the case 
where the growth rate fi(x) is an arbitrary function of the coordinates was investigated by Cantrell 
and Cosner I0. Applying separation of variables to the linearised problem it is possible obtain a 
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system of equations, similar to (13.101 ) 

T' n {t) = X n T n (t) 

(3.17) 

DX'^x) + fi(x)X n (x) = \ n X n (x), 

where X(x) should again satisfy certain boundary conditions. Now, however, fj,(x) ^ const and a 
solution to this equation is known only for some partial forms of the function n(x). If the principal 
eigenvalue of these equations is positive, the stationary solution P(x) = is not stable, implying 
that the population can invade and establish on the habitat. It can be shown that this is always the 
case if the average growth rate is positive, J jj(x)dx > 0. However, this condition is not necessary 
and does not hold even in the simple models considered before. General conditions of uniqueness 
and the existence of positive eigenvalues were obtained by Hess and Kato ||29l , Senn and Hess 
1921 . Cantrell and Cosner g). 

The analysis of 2D and 3D models raises even more questions. In higher dimensions the 
persistence of a population may depend on the geometric form of the favourable patches [15], the 
form of the edges separating the patches and finally it may depend on the behaviour of individuals, 
moving across or along the edges |fT9"1 . 



3.2. Persistence on an infinite habitat 

Travelling fronts While in the previous section we have highlighted some negative aspects of 
diffusion for the fate of a population, in this section we show that depending on the circumstances 
diffusion may as well support population growth. Maybe the most drastic example is the possibil- 
ity to generate the spread or geographic expansion of a population into a new area. Consider, for 
simplicity, an infinite homogeneous habitat which provides a positive growth rate /io everywhere. 
If, as it is commonly assumed (i.e., there is no Allee effect), the stationary state P = is unstable, 
the appearance of organisms in one spot will lead to their expansion over the whole habitat. As- 
suming logistic growth the spatio-temporal dynamics of the population can be represented by the 
following equation 

which is known as the Fisher- Kolmogorov equation after Fisher [20J, who considered the logistic 
dynamics of advantageous genes and Kolmogorov et al. [|48l . who investigated a general form of 
this problem (see also (62l EQ Ml )• 

Fisher ll20l and Kolmogorov et al. [|48l have shown that if initially some part of the habitat is 
not occupied, P(x) = 0, then the population will propagate into this part with the constant velocity 



v f = 2^D. (3.19) 

This solution is illustrated in Fig. [5K which shows a population propagating from left to right 
through a ID habitat. As will be shown below, the magnitude of the propagation velocity is a 
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Figure 5: Front-propagation process in the Fisher-Kolmogorov model (I3.18I) . (A) Propagation of 
a travelling front, if initially P(x, 0) = 0, for x > 0. (B) Propagation of a "pseudo-wave" in the 
limit of zero diffusivity, if initially P(x, 0) = l/(x + l) 5 for x > 0. In both plots P(x, 0) = 1 for 
x < 0. 



crucial factor not only for an invasion process but also for the survival of a population in the 
presence of drift, sinking or other advective processes ||T4l [99l ITOTTl . There is a variety of ways to 
derive relation (13.191) . The more common and rigorous approach suggests to assume a travelling 
solution of the form P(x — vt) and then to prove that this solution is stable only if v = Vf (see e.g. 
ll5"Tl0 . Below we will provide another heuristic derivation confirming the validity of this expression. 

Note that in large aquatic basins the horizontal turbulent mixing increases with the scale of 
phenomena 11731 [671 [6*61 1741 . Petrovskii [[771 1781 showed that this should result in the increase 
of the front propagation velocity. Moreover, this velocity should grow with the size of the area 
occupied by a population. 



Pseudo waves We should stress a condition which is sometimes missed. The velocity Vf is the 
minimal possible propagation velocity, which is realised if the population invades into an empty 
area or if, at least, in this area P(x, 0) decays more rapidly than a Gaussian distribution (see 
Fig. [5K)- Another pattern of so-called "pseudo-waves" may occur for special initial conditions if 
from the start a small amount of biomass is distributed over the whole space. This is illustrated 
in Fig. 03 where the initial distribution of biomass is algebraically decaying for x > and so is 
visually indistinguishable from that in Fig. [5]A As shown this initial configuration yields a much 
faster wave-like spread of the population. Strictly positive initial conditions lead to simultaneous 
logistic growth toward the carrying capacity at every position and, even in the limit of zero dif- 
fusivity, a wave-like pattern arises because the capacity is reached at slightly different times in 
different points and not due to the transport of biomass by diffusion. 

Pseudo-waves appear if the time scale of diffusive transport, td, is slower than the time differ- 
ence Tdem of demographic processes in the neighbouring points 

Ax 2 Ax dP 

r D -^< r dem ~ -j^te » ( 3 - 20 > 

where Tdem can be found expanding the equation P(x, t) = P(x + Ax, t + Tdem) into a series. It 
is easy to check, that for a Gaussian distribution of P(x, 0) both sides of (13.201) are proportional. 
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Figure 6: Schematic presentation illustrating the spread of a population into an oppositely directed 
flow of velocity v. Initially the ID habitat is populated on the left, while the right is not occupied 
by the species. In a system without advection the population would propagate into the right with 
the front velocity Vf, given by equation (I3.19I ). This spread is hindered by the advective flow, which 
is aimed into the opposite direction. As a result the population front propagates with the reduced 
velocity v f = Vf — v. The population can persist only if Vf > v, otherwise it is washed out. 

In contrast, for algebraic or exponential distribution of P(Ax, 0) the right hand side will be larger 
if Ax exceeds a threshold value. Consequently, the diffusive transport will be always slower than 
the time scale r dem of the demographic process, giving rise to a "pseudo-wave". 

Advection Consider now an extension of the Fisher- Kolmogorov equation for a population which 
is additionally subjected to an advective flow with constant velocity v [T63l [641 IX4l [99l ITOTl |4l [50ll 



To investigate the role of advection let us again suppose initial conditions, such that only the 
left part of the habitat is populated, while the remaining part is not occupied by the species (see 
Fig. [6]). In the absence of advection in this system, the population would propagate into the right 
with the velocity Vf. Now assume that this spread occurs in an advective flow, which is going 
into the opposite direction with the velocity v. This can be analysed best by considering a frame 
of reference moving parallel to the flow with the same velocity v, so that in this frame the flow 
velocity is zero. In the moving reference frame the population dynamics obey equation (13.181) 
and the front propagates with the velocity Vf given by (13.191) . Therefore, in the fixed reference 
frame, the propagation velocity is reduced by the velocity v of the advective flow, so that the front 
propagates with the velocity v'* = Vf — v. 

Obviously, the sign of the reduced propagation velocity determines the persistence of the 
species. A negative propagation velocity leads to a population wash-out, whereas a positive prop- 
agation velocity results in the invasion of the empty habitat. Thus, on an infinite habitat in an 




(3.21) 
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advective flow the population can survive only if 

v<v f = 2^D, (3.22) 

or, written in a different form 

v 2 

D > D min = — , (3.23) 

which means that in a flow a population can persist only if the diffusivity exceeds the threshold 
value Dmin- Equations (|3.22l) and (13.231) relate the critical values of diffusivity D and flow velocity 
v and illustrate the constructive interplay that can arise in the presence of both processes, advection 
and diffusion. The region in the (v, D) -coordinate plane, where a population can outgrow an 
advective flow, is visualized in Fig. [7] (area under the solid line). Note that if the habitat is unlimited 
then for any flow velocity v there is a diffusivity which provides the survival of the population. This 
transition plays an important role in many ecological situations and constitutes, for example, a 
necessary condition for the persistence of a population in a river 11991 , as well as for the persistence 
of sinking phytoplankton species in a vertical water column If8~7ll95ll33~l . 



Derivation of the front propagation velocity These results about the spread of a population in 
an advective flow can be elegantly used to derive the population spread, Eq. (13.191 ), in a system 
without advection. For these aims we consider the general model 

dP(x,t) f Nn dP d 2 P /rtnAS 

Here, the local growth term has been linearised for small densities, however we allow for an ar- 
bitrary spatial dependence of the growth rate, ji = fi(x). This model can be simplified with the 
following transformation 

P(M) = exp(||) P(x,t). (3.25) 

The aim of this 'trick' is to eliminate the advective term and as a result we obtain the following 
equation of a system without advection, however with a modified growth term 



dt V 4L> / dx 2 

Note that the transformation (13.251) should be applied also to any boundary conditions of the origi- 
nal system and may alter them. However, if we can neglect the influence of boundaries (e.g. if the 
habitat is infinite, but still heterogeneous) advection obtains a very simple ecological interpretation 
as an additional mortality of strength v 2 /AD. This means that the presence of advection effectively 
reduces the growth rate and this effect increases as D — > 0. 

Now we can "derive" the velocity Vf of the front propagation in the Fisher- Kolmogorov equa- 
tion. On a homogeneous infinite habitat (p{x) = fio) the population can survive only if the growth 
rate is positive. Thus, using (13.261) the velocity should be constrained 

v < 2v^D. 

However, as we showed before, the maximal possible advection velocity, v, is equal to the front 
propagation velocity, Vf, which finally leads us to formula (I3.19I ). 
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3.3. Finite habitats in an advective flow: the "drift paradox" 

While the Fisher- Kolmogorov equation (13.181) assumes an infinite homogeneous habitat, such con- 
ditions are a strong idealisation for most natural populations. To describe some more realistic 
situations, in the following we investigate the role of advection for a population in a heterogeneous 
environment which is additionally constrained by certain boundary conditions. 

Separation of variables Consider again the general model (13.241 ) of the previous section, which 
now is supposed to be complemented by some boundary conditions. Again we use the change 
of variables (13.251) to eliminate the advection term. The linear form of equation (13.261) allows 
a separation of variables P(x,t) = X(x)T(t) and by comparison to Eq. (13.171) we obtain the 
eigenvalue problem for the time-independent eigenfunctions 



while the equation for T(t) remains unchanged. Here, A u denotes the eigenvalues for the problem 
in the presence of advection. Furthermore, X(x) should satisfy certain boundary conditions, which 
depend on the boundary conditions obtained for P(x,t). If we introduce A = X v + v 2 /AD this 
equation will take the form of the second equation of system (13.171) for the problem without flux. 
Thus we can conclude that the presence of advection for a population under boundary conditions 
simply reduces the eigenvalues ||2^l64l[l4]| 



With the same arguments as in Section 13.1.1 to provide the persistence of a population, the 
largest eigenvalue \ v must be positive. Therefore, a population is only able to persist in a flow if 
the proper model without flow has the principal eigenvalue 



KiSS model with advection As a simplest example, consider the KiSS model (13.71) in the pres- 
ence of an advective flow. Recall, that this model describes a population on a finite favourable 
patch provided that the population density vanishes at the boundaries. These boundary conditions 
(P(0, t) = P(L, t) = 0) are not changed by the coordinate transform (13.251) . so that the functions 
P(x, t) and X(x) also must vanish at the boundaries. Therefore, we will obtain the same expres- 
sion (13.1 II) for the eigenvalue spectrum with the dominant eigenvalue A* = /i — n 2 D / L 2 . Taking 
into account (13.281) . we can easily derive the conditions for the persistence of a population on a 
finite favourable patch in an advective flow 




(3.27) 




(3.29) 
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This expression describes a semi-circle in the (D, v) parameter plane (see dashed line in Fig. [7]). 
Note that as L approaches infinity or D approaches zero, the maximal velocity (13.291 ) goes over 
to the condition (13.221) . This means that in the limit of small diffusivity (D <C jiL 2 /ix 2 ) the 
conditions for survival on a finite patch are almost the same as those on a homogeneous infinite 
habitat. Thus, for small D the two critical curves in Fig. [7J almost coincide and show a similar 
behaviour. For larger diffusivity, however, these curves diverge because a finite habitat provides a 
slower propagation velocity (13.291) than an infinite habitat (I3.22I ). Finally, for large values of D the 
increase of the losses across the habitat edges results in the upper diffusivity limit D max on a finite 
habitat. 

Solving this inequality for D, we obtain a limiting interval D min {y) < D < D max {y), with 




1i 2 V 2 



D m in/m,ax ~ — [ 1 ± \ j 1 ~ JY^2 ) ' (3.30) 



'max\ 



Here, D is the critical (maximal) diffusivity (13.131) in the KiSS model. The interval [D min , D n 
specifies the limits of mixing intensity which prevent the population wash-out (D > D min ), but 
still enable the persistence of the population on a finite patch (D < D max ). These values are real 
only if 

v < v max = — . (3.31) 

Note that the critical velocity v max is a threshold when the characteristic time scale of growth 
T gr = 1/pL becomes slower than the time scale of advection r v = L/v. If v > v max the population 
cannot persist, because on the one hand the large advection requires a strong mixing intensity to 
provide the expansion of organisms upstream, but on the other hand such mixing increases the 
transport of organisms into unfavourable areas and the population becomes extinct. 



Critical patch size In the following we investigate the influence of advection on the critical 
patch size in the KiSS model (13.71) and in the model by Ludwig et al. (13.151) . In the first model 
the population vanishes at the habitat edges, whereas in the second the population density should 
vanish when x — > ±oo. Therefore, in both models the transformation (13.251) does not alter the 
boundary conditions and the presence of an advective flow simply reduces the growth rate \i. This 
leads to the effective growth rate 

in the KiSS model and to 

for < x < L 

for x < or x > L , 
in the model by Ludwig et al. 



m 



v 

W 

v 2 
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Figure 7: Parameter regions which permit the persistence of a population in an advective flow. 
Persistence of the population is possible below the critical curves in the (D, u)-parameter plane. 
Two different scenarios are compared. The solid line shows the persistence regime for an infinite 
uniform habitat (cross and diagonal hatching). In this case the limiting velocity scales as the 
square root of diffusivity v ~ \D, equation (13.221) . In comparison, the dashed line shows the 
result for the persistence of a population on a finite favourable patch with lethal boundaries (cross 
hatching). Here, the persistence regime yields a semi-circle in the parameter plane, equation (13.301) . 
Assuming intermediate growth conditions, e.g. similar to the model by Ludwig et al. (13.151) . one 
should expect the critical curve to be located somewhere in-between the solid and the dashed line. 



Substituting the modified growth rate ji" into equation (13.261) . we find similar expressions for 
the critical patch size as in Eqs. (13.121) and (13.161) . where \i should be replaced by fi v . Thus, for the 
KiSS model in an advective flow the critical patch size equals 




(3.32) 



and for the model by Ludwig et al. we obtain 

V = — ^=^= arctan a / — \ — ~*~ , (3.33) 

yjv} ~ V* V 

where Vf = 2^Dfi . Note that both Lq and V approach infinity as v — > Vf. In particular, 
this can happen if D — > D min = v 2 /4/i which is of importance in marine biology as climate 
models predict that the ongoing global warming may result in a higher stratification of the ocean 
water [|6l|9Tl, increasing thereby the requirement on the critical (vertical) patch size for sinking 
phytoplankton species. 
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Persistence in a river Speirs and Gurney ||99ll investigated the conditions for the population 
survival in a river of length L and with a flow velocity v. This problem is known as the "drift 
paradox" because any advection will ensure that the average location of a population will move 
downstream, so at first glance it seems counter-intuitive that a population can persist in a river 
ll28l . In their study Speirs and Gurney used the linear model (13.241) and assumed an impenetrable 
boundary upstream and a totally hostile environment downstream. In this model the critical size L v 
of the favourable patch depends on the ratio between the flux velocity v and the front propagation 
velocity Vf 



Note that L v again increases with an increase of the advection velocity v and approaches infinity 
as v — > Vf. Furthermore, because the favourable patch in this model has only one boundary with 
the hostile environment, we obtain the limit L v — > L /2 as v — > 0, where L is the critical patch 
size (137121) of the KiSS model. 

Extending this model, Pachepsky et al. Il75l derived conditions for the persistence and spread 
of a population of organisms living and reproducing on the sediment and occasionally entering the 
water flow where they can drift and disperse. 

Locally elevated growth rate Dahmen et al. |fl4j considered another extension of the model by 
Ludwig et al. (13.151) . assuming an advective flow and periodic boundary conditions. Furthermore 
they suggested a simple experimental set-up in which a light limited colony of bacteria grows on 
a ring. Almost the whole ring is shaded and unfavourable for the population. Only a small area is 
illuminated through a window of length L, which moves around the ring with a constant velocity 
v. Up to a change of reference frame, this set-up is equivalent to an advective flow and, changing 
the speed v, one can easily regulate the "advection" rate. 

Depending on the parameters and the velocity of the light supply, the population can become 
extinct, can be localised (the maximum of density tracks the location of the favourable patch), or 
delocalised ||64l (due to the periodic boundary conditions the population can survive even if the 
advection velocity v exceeds the propagation velocity Vf provided that the average growth rate is 
positive). Fig. [8] shows an example for such a localisation. 

Solving the eigenvalue problem Dahmen et al. exploited a similarity of equation (13.271) to 
the well studied "square well potential" problem of quantum mechanics ll53l . They showed that 
depending on the dimensionless parameter x = (2/L) y^D/O^o + m), which characterises the 
growth rate in relation to the diffusivity and the habitat length, the critical flow velocity can be 
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Figure 8: Locally elevated growth rate provides the persistence of a population in an advective 
flow with velocity v, Dahmen et al. Ifl4l . A favourable patch = pL > 0) is indicated by the 
horizontal thick line. 
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(3.34) 



where all parameters have the same meaning as in (13.151) . The first solution also corresponds to 
the limit of high mortality (m — > oo) and coincides with equation ( 13.291 ), which was derived for 
the KiSS model with advection. The second solution describes a strongly mixed system. If the 
advection velocity is higher, the population becomes extinct. The value of the critical patch size 
and the critical diffusivity can be easily expressed from equation (13.341) . 

Examining bacterial growth, Lin et al. Il55ll confirmed the results of Dahmen et al. |[T4l ex- 
perimentally and by means of numerical simulation. Joo and Lebowitz ll43l carried out computer 
simulations in a stochastic spatially discrete population model and obtained similar results, con- 
firming the robustness of the model. 



Locally elevated diffusivity Consider a population growing on an infinite favourable patch in an 
advective flow. If the diffusion is too low this patch will not provide a proper propagation velocity 
(13.191 ) and the population will become extinct. Straube and Pikovsky H101II noted that locally 
increased diffusivity can drastically change the situation. A sufficiently large and well mixed patch 
can stop the drift of biomass, stabilising the population dynamics downstream (Fig. [9]). Straube 
and Pikovsky considered the Fisher-Kolmogorov equation with advection (13.211) . assuming that a 
patch of length L has an elevated level of diffusivity D±, which exceeds the diffusivity D in the rest 
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Figure 9: Locally increased diffusivity stops the population wash-out in an advective flow with 
velocity v, according to Straube and Pikovsky H101II . The horizontal thick black line indicates the 
location of the patch with increased diffusivity D 1 . 



of the habitat. Note that if the diffusivity depends on the coordinate, D = D(x), equation (13.211 ) 
takes the form 
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Straube and Pikovsky H101H showed that the population survives, if the size of the intensively mixed 
patch is larger than 
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As the diffusivity D\ approaches infinity the critical patch size approaches 



lim L r 



v 2 — v 2 



Thus, even for the infinite mixing intensity the critical patch should be of finite size. 



Techniques from quantum mechanics Birch et al. [4] considered the Fisher-Kolmogorov equa- 
tion with variable growth rate and advection on a 2D plane. In this study they made use of the struc- 
tural similarity between equation (|3.27l) and the time-independent Schrodinger equation. Based on 
this similarity, they demonstrated a few examples where the application of perturbation theory, the 
method of Wentzel, Kramers, and Brillouin (WKB), and other quantum mechanical techniques are 
beneficial for the analysis of equations similar to (I3.27I ). In particular, Birch et al. were able to de- 
termine a trade-off between the critical diffusivity and growth rate, which provide the persistence 
of the population in this system. We believe that the application of such well-known techniques 
from quantum mechanics has not been fully exhausted yet. Such methods could be a promising 
direction for the further development of the theory of extended populations. 
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Summary At the end of this section we want to summarize the main results. On an infinite 
homogeneous patch a stationary state is achieved when the growth rate is equal to the mortality. 
On a finite patch, which adjoins some unfavourable environment, the growth rate should exceed 
mortality to compensate for losses across the patch edges. The internal patch production grows 
with the patch size and with the growth rate, whereas the losses at the edges increase with the dif- 
fusivity. The spread of organisms over a habitat is defined by the propagation velocity of the front 
of the population density. This velocity scales as a square root of the growth rate and diffusivity. A 
population can survive in a flow only if the propagation velocity is higher than the velocity of the 
flow. 

Thus, an increase of the growth rate improves the conditions for the survival on a finite favourable 
patch and in a flow. The increase of the patch size increases the internal patch production, whereas 
the losses across the patch edges show only a weak dependence. The critical patch size is achieved, 
when the internal production is balanced by the external losses. The size of a large habitat has 
a small influence of the propagation velocity. However, an unfavourable environment around a 
small favourable patch truncates the propagating front, this can essentially reduce the propagation 
velocity and worsen the conditions in a flow. Finally, an increase of diffusivity on a finite habitat 
increases the losses into the unfavourable environment, which can lead to the population extinc- 
tion. However, a minimal level of diffusivity is necessary to prevent the population wash-out in a 
flow. These two opposing processes result in a diffusivity window, which provides the population 
survival on a finite habitat. This window exists only if the flow velocity is less than some critical 
value. 

4. Vertical phytoplankton distribution 

In the previous chapter, using simple one-dimensional models, we discussed the influence of mix- 
ing and advection on the population survival in a heterogeneous environment. To model a non- 
uniform environment, we simply used an explicit form for the growth rate, fi(x), however we did 
not discuss the origin of this heterogeneity. In this chapter we extend these results by considering 
more complex models which describe resource limited population growth. These models are based 
on simple physical and biological laws and consistently describe the dynamics of a population and 
its limiting resources. Thus, they demonstrate natural mechanisms which lead to the appearance of 
favourable patches. The intention of this chapter is twofold: on the one hand we aim to illustrate 
how the main findings of the previous chapter apply to the context of consumer-resource models. 
On the other hand we show new effects arising due to new properties, which are not present in the 
simple models. 

The main object, which we will use for illustration, is the dynamics of a vertical phytoplankton 
distribution. This is important as phytoplankton are the primary producers in almost all aquatic 
food webs with a major influence on nearly all freshwater and marine ecosystems. The two main 
factors limiting the production of phytoplankton are the availability of nutrients and light. To un- 
derstand how these resources affect the phytoplankton biomass, consider their distributions in a 
water column. In general, the light intensity reduces with depth and, in nutrient rich regions of the 
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ocean, the well illuminated surface layer constitutes a favourable area for photosynthetic phyto- 
plankton species. By contrast, the nutrient concentration can behave just in the opposing way. The 
sedimentation of dead biomass (detritus), with the successive remineralisation in the deep layers or 
in the sediment lfT6l causes an increase of the nutrient concentration with depth H117II . Thus, while 
light limitation may lead to the formation of a surface phytoplankton maximum, a lower nutrient 
concentration favours a phytoplankton build-up in deeper layers. This tension between light and 
nutrient limitation from two opposite sides frequently causes optimal growth conditions in subsur- 
face layers (e.g. 03 \T3[ [3T]|). This fact often leads to the appearance of maxima of chlorophyll 
or biomass distributions at approximately 30-100 m depth. So called deep chlorophyll maxima 
(DCM) (D [HH HUB ED and deep biomass maxima (DBM) flUld are ubiquitous phenomena 
and can be observed in many oligotrophic regions in the ocean, marine systems, and deep lakes. 

Another important component of a stratified water column is an upper mixed layer (UML). A 
UML commonly occurs in oceans and lakes due to mechanical perturbation of the surface waters 
(e.g. due to wind, waves, and storms). This layer is separated from the deep layers by a thermocline 
[fT2l . which is defined as a relatively thin layer below a UML characterised by an strong change 
in temperature with depth. Mixing in a UML is much stronger than in the layers below it. As a 
result, the distributions of nutrients, temperature, salinity, etc. are nearly uniform in a UML and 
have gradients below it. The depths of a UML can usually vary from 10 m to 100 m, see e.g. H110H . 

Compared to the models of the previous chapter, the system behaviour in a water column can 
be further complicated due to a feedback loop between the biomass and resource distributions. The 
growing biomass shades light, consumes nutrients and is remineralised, which ultimately changes 
the total resource distribution. This, in turn, can lead to a new biomass distribution, which will 
generate a new resource profile and so on. As will be shown below, these complicated, self- 
organised dynamics can lead to new phenomena and diverse behaviour. For example, if the mixing 
is small, the final solution becomes non- stationary and oscillates 11361 . whereas in the presence of 
an upper mixed layer the system may exhibit bistability and the solution may be sensitive to the 
initial conditions H118ll89ll . 



Equation of growth To formulate a mathematical framework for this chapter, let us consider a 
vertical water column of depth Z B . Let P(z,t) denote the density of phytoplankton at time t and 
depth z. Note that z = denotes the sea level surface and the z-axis is directed downward. For 
the sake of simplicity, assume that phytoplankton growth is limited only by the availability of light 
and a nutrient (the model can easily be extended to take into account multiple nutrient limitation 
(261). hi our approximation the dynamics of a phytoplankton population obey a reaction-diffusion- 
advection equation, similar to equation (13.11 ) considered in the previous chapter (see If85ll95ll46ll36l 
among others) 

where m is the mortality (compare to equation (12.31) ). v is the phytoplankton sinking velocity, and 
D is the diffusivity, which in general can depend on z. 

Furthermore, the growth rate fi(N, I) depends on the local values of light intensity I(z, t) and 
nutrient concentration N(z, t) at each vertical position. If both nutrients are essential, n(N, 7) can 
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be represented in the form of Liebig's law of minima 

^N,I) = f i mm[f I (I),f N (N)} . (4.2) 
It can be also written in the multiplicative form 

li(N,I) = nofi(I)f N (N). (4.3) 

Here fj, is the maximum growth rate and //(/) and /at(-/V) describe the limitation by light and 
the nutrient. The specific form of these functions depends on many factors, for instance, strong 
light may photoinhibit photosynthesis and reduce the growth rate for large values of I If83ll22~ll86ll . 
However, usually it is suggested that f(x) — > 1 as x — > oo, that is, the maximum growth rate is 
achieved when all resources are unlimited. For phytoplankton modelling, the most frequently used 
form is the Monod (or Michaelis-Menten) kinetics H107H 

fl{1) = , MN) = , (4.4) 

where Hj and H N are the half-saturation constants for nutrient-limited and light-limited growth, 
respectively. However, this non-linear form often admits only numerical investigation. Analytical 
solutions are commonly possible only for a linear or algebraic form of f(x) [|95l[T7l . 



Boundary conditions By default, we assume that the surface and bottom are impenetrable for 
phytoplankton 

' dP 
vP(z,t)-D 



dz 



= . (4.5) 

o,z B 



To model a stratified water column, one can either separately solve the equations in a UML 
and below it, supposing infinite mixing within the UML and a small diffusivity Do in deep layers. 
Assuming continuity of the flux across the thermocline, we obtain the boundary condition at the 
bottom of a UML (see e.g. 0511 ) 



dP\ 

\ z =z T -o= [vP(z)-D D ^\ 



=Z T +0 



where Z T is the depth of the thermocline. On the other hand, to simulate the water column in a 
single framework, one can assume a gradual transition from a UML to the deep layers ll89ll 



D & = D ° + i + e(*-f T V« ' (4 ' 6) 
where Djj and Dd are the diffusivities within and below a UML, respectively, and the parameter 
w characterizes the width of the thermocline. 

So far, we did not specify any equations for the distribution of nutrients N(z,t) and light 
I(z, t). In the next two sections we will review several models which couple the light and nutrient 
dynamics with phytoplankton growth (14.11) . First we will consider models in which phytoplankton 
growth is limited only by the light availability, whereas in a second step, we will review more 
complex models, incorporating both light and nutrient limitation. 
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4.1. Light limitation 

In this section we will consider theoretical models in which the light gradient is a key factor. Such 
models give an adequate description for eutrophic aquatic environments (as observed in many re- 
gions), where the nutrients are in ample supply and light becomes a crucial factor which determines 
the distribution and the dynamics of phytoplankton Il59ll90l [TTTl . So in the following we assume 
that the nutrient dependence of the growth rate is saturated, /n(N) — > 1, and we can neglect the 
limitation of growth by nutrients. 

The spatial profile of light intensity in a water column is described by Lambert-Beer's law (see 
e.g. 11451 ) which states that the gradient of light intensity at depth z is proportional to the light 
intensity at this depth 

f = -*I . (4.7) 

az 

The coefficient k includes both the absorption of light by water and the attenuation by the phyto- 
plankton cells 

k = K bg + kP(z,t) , 

where K bg is the background turbidity and k is the per capita attenuation coefficient of the algae 
cells. Integrating (14.71) from surface to depth z, we obtain 



I(z) = /i„exp 



-K bg z- [ kP(t,z')dz' 
Jo 



(4.8) 



where I in is the light intensity at surface. 

Equations ( 14.11 ) and ( 14.81) , being coupled by means of the growth rate ( 14.21) or ( 14.31 ), yield an 
integro-differential system of equations. It is not straightforward to obtain rigorous or analytical 
results for such a system and even a numerical solution encounters certain difficulties If35l 1 10311 . 
Nevertheless, without solving any equations, it is clear that the light intensity in the water column 
is reduced with increasing depth. Thus the light limitation forms a favourable area close to the 
surface, and the dynamics of the phytoplankton population should be related with the results of 
the previous chapter, obtained for heterogeneous environments in the presence of advection and 
diffusion. 



Critical values for phytoplankton growth Depending on the depth of a water column or on the 
diffusivity, a light limited phytoplankton population can survive or become extinct. Huisman et al. 
113311351 combined the conditions for survival into a single conception of the critical conditions for 
phytoplankton blooming in a closed water column (Fig.fTOK) and within an UML (Fig.[T0B). The 
main difference between these models is that in a closed system, the sinking of cells is stopped at 
the bottom, whereas in an UML the biomass can sink across the thermocline to the deep aphotic 
layers. 

To determine these conditions, consider first an unstratified water column with a constant dif- 
fusivity and assume impenetrable boundary conditions ( 14.51) for the phytoplankton biomass. Using 
the system of equations (14.11 ) and ( 14.81 ) in the limit of zero background turbidity, K bg = 0, and for 
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Figure 10: Critical conditions for phytoplankton blooming: (A) in a closed water column and (B) 
in an upper mixed layer. Both figures are reprinted from Journal of Sea Research, 48, Huisman 
and Sommeijer, pp. 83-96 11351 . Fig. 4 and 6, with permission of the author. 

a general monotonic growth rate Shigesada and Okubo ll9~5l showed that a sinking phyto- 
plankton species can establish a population only if 

This expression coincides with condition (13.231) from the previous chapter and implies that the 
minimal diffusivity should provide a front propagation velocity which is larger than the sinking 
velocity. The same condition was derived earlier by Riley If8"71 and other authors from the Fisher- 
Kolmogorov equation. It is interesting to note 11951 that if K bg = and a non-trivial solution 
exists, then the total biomass in this model does not depend on the sinking velocity. The sinking 
just shifts the bulk of biomass downward, preserving, however, the total amount of biomass in the 
water column (Fig.fTTK). 

Ishii and Takagi 11401 relaxed the condition Kb g = and proved some existence, stability and 
uniqueness results for this system. Assuming an algebraic form of the growth rate, ~ I a , 
Ebert et al. IfTTl have found some approximation for the minimal diffusivity, D min , and for other 
critical parameters. 

If a water column is sufficiently deep and K^g ^ then the net production rate is positive only 
above the compensation depth, Z c , which is defined as the depth z at which the local production 
rate is zero in the absence of biomass (Fig.fTTB). From the light attenuation curve (14.81 ) we find 

In I in - In I c 

— — 



K, 



where I c is defined as the compensation light intensity at which the growth rate is equal to mortality, 
= m, thereby the compensation depth is species specific. 
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Figure 11: (A) Vertical profiles of phytoplankton for different values of the dimensionless sinking 
velocity u = v/fj, D ||95l . While the shape of the profile is changed, the total amount of biomass 
in the system remains unchanged. The figure is reprinted from Journal of Mathematical Biology, 
12, Shigesada and Okubo, pp. 31 1-326 Il95l , Fig. 4, with permission of the author. (B) Schematic 
representation of the compensation depth Z c and the critical depth Z cr . 



If a water column is shallower than the compensation depth and we assume that the bottom is 
impenetrable for the biomass, then the population can survive even if mixing is less than a minimal 
diffusivity (14.91) because the cell settling will be stopped at the bottom (Fig.PTOR). 

In general, the compensation depth divides a water column into a favourable and an un- 
favourable regime (Fig. [TTB). In a well mixed water column the losses in the deep layers can 
lead to the population extinction. However, they can be compensated by the production in the 
euphotic zone, if the unfavourable region is relatively small. Considering a simple mathematical 
model of a well mixed water column Sverdrup HI 0211 defined a critical depth as the depth of a water 
column at which the total growth is equal to the total loss of biomass. Similar to the compensation 
depth, the critical depth can be reinterpreted in terms of the critical light intensity ll37l 

In I in - In I out 

K bg 

where I out is the light intensity at the bottom of a sufficiently shallow (Z B < Z cr ) closed water 
column after the light limited population of phytoplankton have reached an equilibrium state. The 
critical depth depends on many parameters, it increases with the incident light intensity and with 
the phytoplankton growth rate, and it decreases with the mortality rate Il3~7l . 

In a well mixed water column, an excess of the critical depth over the compensation depth 
determines the maximal possible losses in dark layers, which can be still compensated by the 
production in the euphotic zone. However, similar to the KiSS model or to the model by Ludwig 
et al. (see Sec. 13.1.1) these losses diminish with a decrease of mixing. Extending this research, 
Huisman et al. 11371 showed that if the depth of a water column or a thermocline exceeds the critical 
depth, the population survival still is possible if turbulent mixing is less than a maximal diffusivity. 
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This critical condition is similar to the existence of the maximal diffusivity (13.131) and (13.301 ) in the 
KiSS model, however here it describes the behaviour of a more realistic system. Therefore similar 
to the KiSS model in advection (13.301) , a sinking population can survive in a water column of any 
depth if mixing remains between a minimal and a maximal value (Fig.[l0lV). 

We now turn to a stratified water column with a UML. We assume that the mixing in the 
deep layers is less than a critical value (14.91 ), so that the population survival will depend on the 
characteristics of a UML. Condie and Bormans |[T2ll showed that if a UML is shallower than 

„ v 



a population cannot survive (compare with (13.311) ). In other words, for the survival in a UML, 
the demographic time scale should be faster than the characteristic time of advection. However 
usually, Z T ^ min is sufficiently small and this criterion is satisfied. In this case the population can 
persist if the strength of mixing remains within the turbulent window [D min , D max ] (see Fig.flOB). 
Furthermore, if the diffusivity exceeds a maximal value the population survives if the depth of a 
thermocline is smaller than a maximal depth, which is defined as the maximal depth of a well mixed 
upper layer at which losses and production are equal. This depth is slightly smaller than the critical 
depth in a closed water column, owing to additional losses of biomass across the thermocline. 

The fact that a deep upper layer can prevent phytoplankton blooming was noted experimentally 
in 1935 by Gran and Braarud E3l . who investigated the conditions of phytoplankton blooming in 
the upper mixed layer. They reported that until there exists a deep UML, phytoplankton production 
cannot exceed the destruction by respiration and phytoplankton blooming is not possible. The 
concept of the maximal diffusivity is also consistent with field experiments, see e.g. HI 061 [18117611 . 



4.2. Light and Nutrient limitation 

In the last section we will discuss models which take into account both light and nutrient limita- 
tion of phytoplankton growth. These models are more difficult to analyse and often admit only 
numerical investigation. However, they are more realistic and provide some understanding of the 
processes occurring in deep waters of many regions where surface layers are nutrient depleted 
If831lll3ll82llll51lll21ll09ll36ll . Furthermore, in the tension of two opposing resource gradients 
the location and the size of a production layer becomes a function of the phytoplankton abundance 
and the initial conditions, that can lead to new patterns and new dynamical behaviour. 

A coupled system of reaction-diffusion equations describing nutrient-phytoplankton cycling 
was probably first investigated by Okubo ll69l ITU . Radach and Maier-Reimer [[851 suggested 
a mathematical model of phytoplankton growth which included light-nutrient-phytoplankton dy- 
namics. This model was extended by Jamart et al. pT| who considered limitation by two nutri- 
ents, grazing and the variability of the parameters with depth and time. This approach (see also 
111131 l82l 11151 111210 gave rise to a growing set of ecological models, which include cycling of 
many chemicals H117H . coupling with meteorological data 11421 . interplay of different phytoplank- 
ton groups, and 3D simulations ff60~ll2TTl . Here, however, we will focus on the theoretical aspects 
and consider only simple conceptual models. 
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Conservative models The nutrient dynamics include uptake by phytoplankton, remineralisation 
of dead biomass back into a nutrient pool and diffusion. Assuming absolutely effective recycling 
we obtain 

- 

(4.10) 

dN(z,t) , %T r s n „ ^d 2 N 

sit = -WW + + 

where the biomass is measured in terms of its nutrient content (compare to the non-spatial version 
(12.31) of this model). We do not include advection in the second equation, as nutrients, which are 
dissolved in water, are only slightly influenced by the gravity force. Nevertheless, this term should 
appear, if advection is caused by a vertical or horizontal stream. 

Furthermore, we assume that the nutrient cannot diffuse across the surface and a large nutrient 
pool in the sediment or in deep ocean layers sustains a constant concentration, Nb, the bottom of 
the water column 

™jjh$- = 0, N(Z B ,t) = N B . (4.11) 

Fig. E2] shows typical final distribution of phytoplankton and nutrient given by model (14.101) . 
supplemented by equation (14.81) for light. Hodges and Rudnick [|30l pointed out that, independent 
of the functional form of the growth rate and of the light distribution (assuming that light decreases 
with depth), this model can reproduce a deep stationary phytoplankton maximum only if v > 0. 
In other words, the presence of opposing resource gradients is not sufficient to reproduce a deep 
phytoplankton maximum. To prove this, let us define the total concentration of the nutrient as 
S = P + N. Consider an equilibrium state, when the left-hand-side of (14.101) equals zero. By 
adding both equations (14.101) we obtain 

D &S_ v dP = Q 

dz 2 dz 

Assuming v — and integrating this equation over z we find 

OS OP ON n 

— = const = — + — = , 

oz oz oz 

owing to the boundary condition (14.51) and (14.1 II) at the surface. Thus S = N + P = const and a 
deep phytoplankton maximum should be accompanied by a deep minimum of nutrient. However, 
if the light intensity reduces with depth, this profile is unstable because there is no factor limiting 
phytoplankton growth in the upper layer. Thus this system should exhibit a surface maximum 
(Fig.[T2"R). However, similar to the model without nutrient limitation (Fig.fTTK). the phytoplankton 
sinking shifts the maximum of biomass downwards (Fig.[T2B). 

Extending this model, Hodges and Rudnick |30| included a detrital pool, T(z, t), as the third 
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Figure 12: (Colour online) Typical distributions of phytoplankton (green), nutrient (blue), and 
light (red) in the conservative model (14.101) without sinking (A) and with sinking (B), according to 
Hodges and Rudnick [|30l . 
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(4.12) 



where vp and vt are the sinking velocities of phytoplankton and detritus respectively. Note that 
usually detritus sinks much faster than phytoplankton |[94l[84l . In this model the cycle of chemicals 
includes three stages: the transfer of biomass to detritus with mortality m, the remineralisation of 
detritus back into nutrients with remineralisation rate r, and finally the consumption of nutrient 
by biomass. While this model can exhibit deep maxima if v = 0, the change of phytoplankton 
concentration is very small and cannot represent real data. An apparent maximum can be observed 
only if one assumes sinking of detritus or phytoplankton. Hodges and Rudnick extended this state- 
ment to any number of compartments, which however do not include depth dependent parameters. 
Thus, sinking is a major component of this system. The sedimentation of organic matter removes 
the nutrient fixed in phytoplankton cells from the upper layer, which leads to the formation of deep 
phytoplankton maxima. 

Beckmann and Hense [3 ] performed numerical simulations and analytical evaluations of model 
( 14.121) , assuming that detritus sinks relatively fast, whereas the phytoplankton sinking is negligible. 
Fig. |T3l reproduces a typical distribution of physical characteristics in this model. Furthermore, 
Beckmann and Hense suggested to extend the concept of compensation depths. Instead of the static 
definition in the absence of biomass they suggested to use two dynamical depths at which the in situ 
production rate of phytoplankton is zero, owing to the light or nutrient limitation. If phytoplankton 
sinking velocity is zero then in equilibrium (see (14.121) ) these values can be expressed from the 
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Figure 13: (Colour online) Typical distributions in the three compartments model (14.121) . (A) 
Vertical profiles of phytoplankton (P), detritus (D), nutrient (N), and light (/). (B) The upper 
and lower limits of the production layer (dashed lines) are the species compensation depths. The 
figures are reprinted from Progress in Oceanography, 75, Beckmann and Hense, pp. 771-796 
Fig. 2, with permission of the author. 
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where zi N ^ and Zc > are the compensation depths due to nutrient and light limitation (Fig.PT3~B). 



Non-conservative models Model ( 14.121) contains three reaction-diffusion equations and an equa- 
tion for the light distribution. This makes further analysis difficult. However, since detritus sinks 
relatively fast [|94ll84l . we can simplify the model assuming that a part e of the dead biomass is 
instantly remineralised in situ, whereas the rest sinks until it reaches the bottom and sustains a con- 
stant nutrient concentration at the bottom. Thus we obtain the following non-conservative system 
of equations 



dP{z,t) 
~~dt 

dN(z,t) 
di 



dP d 2 P 
H(N, I)P-mP-v— + D— 



(4.13) 
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d 2 N 
dz 2 
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This model (compare to equation (12.71) ) can reproduce deep phytoplankton maxima even if the 
sinking velocity is zero, owing to the fact that a part (1 — e) of the fixed nutrient is implicitly 
transferred from the upper layer to the bottom. Even though this model is non-conservative and 
has apparent disadvantages, it or similar models were successively applied to reproduce field data 
iflTl 11091 l36l II 181 . However, we are not aware of any comparison of the two model classes (14.101) 
and (14.131) . which might be interesting. 

In the case of zero sinking velocity, Klausmeier and Litchman 11461 performed analytical calcu- 
lations for model (14.131) . Assuming that the phytoplankton distribution can be approximated by a 
Dirac 5-function and further that an infinitely small production layer should be located to balance 
the light and the nutrient limitation, Klausmeier and Litchman found an equation for the position 
Z* of a deep maximum, which for boundary conditions (14.51) and (14.1 II) reads as 

In (J m /J c ) K bg HqD{N b - N c ) 

k k m(l-e)(Z B - Z*) ' 

where N c and I c arc the critical values of light and nutrient intensity for which the growth rate is 
equal to the mortality rate. 

Oscillations and chaos Huisman et al. Il36l pointed out that system (14.131) exhibits oscillations 
of biomass if the mixing is reduced below a critical value. Fig.[l4]shows the behaviour of biomass 
and of nutrient in two typical cases. In the first case (Fig. [14a) the mixing intensity is high enough 
to provide a stable distribution of biomass. If however the level of diffusivity is reduced, then only 
oscillatory, or even chaotic patterns, can appear (see Fig. [T4~b and [T4~b). As noted by Huisman 
et al. these oscillations are caused by the difference in the time scales of the rapid transport of 
phytoplankton, consuming the nutrients, and the slow upward transport of nutrients. Furthermore, 
as shown in Section l33l for the survival of a population in an advective flux the diffusivity should 
exceed a minimal level (13.301) . which increases with the reduction of the habitat and of the growth 
rate. In the absence of biomass the nutrient can be nearly uniformly distributed over the water 
column, thereby the growth rate becomes only light limited and the production layer extends from 
the surface to the compensation depth, which is usually sufficiently large. Thus, without biomass, 
the level of mixing might be sufficient to induce population growth. However, the consumption 
of nutrients and self- shading of light reduce both the growth rate and the width of the production 
layer. That, in turn, increases the value of the minimal diffusivity (13.301) and finally the sinking 
may lead to the population wash-out if the diffusivity in the water column becomes insufficient. 

Koszalka et al. Il50l noted that the periodic oscillations of phytoplankton biomass will be most 
probably disguised by currents and horizontal inhomogeneity in a real ecosystem. 

Upper mixed layer Hodges and Rudnik ll30l and Beckmann and Hense showed that if self- 
shading of light can be neglected in equation (14.81) . then an upper mixed layer does not lead to any 
qualitative changes in the system dynamics. However, Yoshiyama and Nakajima fl 1 1 811 pointed out 
that a UML can lead to bistability of phytoplankton profiles. 

Ryabov et al. Il89l generalised this result by taking into account the competition of two species 
and relaxing other assumptions. They considered the model (14.131 ). assuming a gradual change 
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Figure 14: (Colour on line) Evolution of the phytoplankton density and the nutrient concentration 
with time, (a) A stable DCM (D = 0.5 cm 2 /s), (b) small oscillations in the DCM (D = 0.2 cm 2 /s), 
and (c) large-amplitude oscillations in the DCM, with double periodicity (D = 0.12 cm 2 /s). The 
figures are reprinted from Nature, 439, Huisman et al., pp. 322-325 ll3"6ll . Fig. 2, with permission 
of the author. 

of diffusivity <|4.6I) from a UML to the deep layers, and showed that under certain parameters, 
depending on the initial conditions the production layer can be steadily located either within a 
UML or below it. Fig. [15] provides a rough insight into the system dynamics. In the absence of 
an upper mixed layer the difference in the locations of biomass and of the production layer drives 
the bulk of biomass towards the production layer (Fig. [T5K). The shift of biomass can lead to the 
redistribution of resources, which in turn can change the location of the production layer. This 
process repeats until the system reaches an equilibrium configuration (Fig. [T5B), when the centre 
of biomass coincides with the centre of production. Now consider a system with an UML. In a 
certain range of parameter the UML does not affect distributions with a deep maximum of biomass 
(Fig.[T5C). However, the initial growth of biomass within the UML begets another stable solution 
with a maximum of biomass located within the UML. The biomass is almost uniformly distributed 
within the UML and its location is uncoupled from the location of the production layer (Fig. [150). 
As a result, a gradual shift of the bulk of biomass into deep layers is no longer possible and the 
transition to a deep biomass maximum can only take place if the light intensity below the UML 
is sufficiently large to provide positive net growth in deep layers - otherwise the phytoplankton 
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Figure 15: Typical vertical phytoplankton profiles P(z) in a system without a UML (top) and with 
a UML (bottom), assuming a gradual change of diffusivity (14.61) in model (14.131) . Without a UML, 
a non-stable phytoplankton distribution (A) evolves to a single stable solution (B). Under the same 
conditions in the system with a UML, we observe two stable distributions: with a maximum in the 
deep layers (C) and with a maximum in the UML (D). The dot-dashed and dashed lines show the 
limitations of growth (14.41) by light and by nutrient, respectively, vertical dotted line shows the level 
of mortality. Black and grey arrows show the centers of biomass and net production, respectively. 

remains trapped in the UML. Thus the production layer can occupy different parts of the water 
column, depending on the current system state and on initial conditions. 

5. Discussion 

Concluding this review we would like to make a few notes. First, let us compare the behaviour of 
the critical patch models and those based on the consumer-resource dynamics. The latter models 
can be divided into two large groups. In the first group we would include those systems in which 
the location of a favourable patch is constrained by some environmental conditions. For instance, 
the limitation of phytoplankton growth by light leads to the formation of the favourable patch in 
the upper level of a water column. The dynamics of this group and of the critical patch models 
demonstrate many general traits and many effects can by predicted and evaluated on the basis of the 
minimal models. The second group consists of the models in which the location of the favourable 
patch is determined by the dynamical interplay of different factors. For example, we can consider 
the growth of phytoplankton biomass driven by two opposing resource gradients. In this group, 
the location of the favourable patch is not predefined. Moreover, the system dynamics becomes 
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very sensitive to the implementation of the consumer-resource cycling. This complexity leads to 
the arising of new patterns and new dynamical behaviour, which can hardly be reproduced in the 
framework of the critical patch models. 

The second remark concerns the advantages and disadvantages of partial differential equations 
(PDEs) for modelling ecological systems. PDEs provide very convenient and powerful tools for the 
investigation of population dynamics. First, in the same framework, we can consider such different 
and complex phenomena as, for instance, the vertical distribution of sinking phytoplankton cells or 
the survival of a population drifting in a flow. Second, analytical solutions in many cases provide 
important predictions and understanding of the main effects, which can appear in more realistic 
systems. Third, one can perform an exhaustive numerical simulation of a model, determining all 
possible bifurcation points. Finally, seemingly the pool of methods developed for the analysis of 
partial differential equations is not played out yet and this approach can still gain a lot of useful 
techniques from quantum mechanics and statistical physics. However, we would like to mention 
as well some restrictions of this approach. Intrinsically it is always suggested that this approach is 
suitable for systems containing many organisms, so that the relative fluctuations of density become 
negligible and all function are continuous. However this statement does not hold if we consider the 
survival-extinction transition. As the system approaches its critical state, the population density 
declines and the fluctuations of density (demographic stochasticity) start to play a crucial role 
H1051I . Thus, in reality, the extinction of a population might occur under conditions which still 
allow for the population survival in a deterministic PDE framework. Therefore, the development 
of a theory including stochastic effects is necessary for the correct representation of the transient 
behaviour. 
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